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Abstract 

During the past three years, Wapenaar, Snieder, Broggini and others have developed an 
algorithm to compute the Green's function for any point inside a medium to points on the 
surface from measurements on that surface only. Their algorithm is based on focusing an 
incoming wavefield to a single point in order to create a virtual source at the focus. The 
procedure has been justified only by heuristic arguments. In this paper I am using simple 
physical arguments to prove an integral equation for single-sided, higher-dimensional inverse 
scattering. This integral equation is equivalent to the Wapenaar iteration algorithm. The 
equation will be exact, including all internal multiple reflections. The derivation makes use 
of time-invariance but does not use the explicit form of the wave equation. It is therefore 
not only applicable to the acoustic wave equation, but also to other time-reversal invariant 
systems. Potential areas include seismology, electronics, microwave and ultrasonic inverse 
scattering as well as quantum physics. The simplicity and generality of the argument will 
make this paper accessible to researchers from all the mentioned fields. 



1. Introduction 



Recently Waapenar, Snieder, Broggini and others dug out work by Rose p!P|5] on one- 
dimensional inverse scattering, extended it somewhat to be able to read off the Green's 
function and generalized it to three dimensions [5|6f7f8] . The idea is as follows. Given sources 
and receivers on the surface of a medium, a shot gather is recorded. Using the data of the shot 
gather, one iteratively constructs an incoming wavefield emitted at the sources, which focuses 
at a time t = entirely at a single point x in the subsurface. 'Focusing' means that at a 
specific point in time, the wavefield vanishes everywhere except at x. Since the basic acoustic 
wave equation is invariant under time-reversal, knowing the focusing wavefield for some point 
X amounts to knowing a wave originating at that point and propagating backwards towards 
the sources. This is the Green's function of that point, called the virtual source. The appeal 
is that apart from noise and numerical inaccuracies, the method is exact and takes into 
account all internal reflections. A velocity model is only needed to determine the position of 
the virtual source but not as an input to the calculation of the Green's function. The Green's 
function for every point in an acoustic medium encodes its entire information and can for 
example also used for seismic imaging as was done in [10] . The iterative algorithm for ID had 
been developed in [Ij and amounts to solving the Gelfand-Levitan equation - also called 
Marchenko equation or Gelfand-Levitan-Marchenko equation. In [5J the theory was applied 
to ID synthetic data and in |7j the aforementioned iteration algorithm has been generalized 
to three dimensions and justified by heuristic arguments. Having a more rigorously proven 
three-dimensional analog to the Gelfand-Levitan equation would be beneficial for a number 
of potential applications in areas such as seismic interferomety, seismic imaging, inverse 
scattering of microwaves, of ultrasonic waves and possibly scattering in quantum mechanics. 
The purpose of this paper is to obtain the exact equation for single-sided inverse scattering 
from simple physical arguments and very basic math. The resulting equation is equivalent to 
the aforementioned iteration scheme. Furthermore, since inverse scattering equations have 
been used in a variety of areas and contexts, a literature search is at times confusing. I intend 
to clear some of the fog by gathering the known scattering equations and their application 
in the last section of this paper. 

2. Derivation of the inverse scattering equation 

The reflection of a wave incident on a medium is a wavefield R{x, t) propagating in the 
opposite direction to the incident wave. In practice, the wave is recorded at a receiver at 
a fixed position and we assume that beyond the receiver it keeps propagating freely so we 
can write R{x + t) just from measurements at the receiver. In one dimension, we assume 
the source and receiver to coincide. In higher dimensions, a source s G anywhere on the 
medium surface S causes a scattered signal to arrive at all receivers g G 5*. Therefore the 
reflection response in a higher- dimensional setting has one more parameter, R{g,s,t). In 
geophysics this is the shot gather which is recorded in seismic experiments. 
Physically speaking, the Gelfand-Levitan equation connects the reflection data from the 
scattering problem to the wavefield of a focusing wave. Numerous derivations of the ID 
Gelfand-Levitan equation for inverse scattering already exist, for example [2|3P] . Here the 
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(a) Wavefield for the reflection problem: 



X ^ M 


X e M 


5{x -t) + R{x + t) 


K{x, t) (with K{x, t) = for t < x.) 


(b) Wavefleld of a focusing wave: 




X 


X e M 


6{x-t) 


5{x-t)+ g'{x,t) 


(c) Wavefleld of the reflection problem in terms of focusing waves: 


X ^ M 


X e M 


1 6{x - t') [6{t - f) + R{t + t')] dt' 

J — oo 


j 5{x - t') [5{t - t') + R{t + t')] dt' 


=5{x-t)+R{x+t) 


+ 1 g'{x, t') [5{t - t') + R{t + t')] dt' 

J — oo 



Table 1: In (c) we superimpose focusing waves in such a way that the field is identical to (a) 
in the region x ^ M. 

equation is rederived in a way which allows generalization to higher dimension. The basic 
idea is simple: We write down the wavefields of both the refiection problem as well as the 
focusing wave. Since the wave equation is linear, we can superimpose solutions, including 
solutions shifted in time. We superimpose focusing waves in such a way that the resulting 
wavefield is brought in agreement with the scattering problem. The result will be the desired 
integral expression which we can solve for the unknowns. That's all there is to it. 

2.1. Derivation in ID 

In our setup of the ID scattering problem, an impulse signal 5{x — t) is incident from the left 
on a medium M. One part of the wave is refiected, R{x+t), while another, as of yet unknown 
part K{x,t), enters the medium. Typically the refiected wave can be measured while the 
field K{x, t) is unknown. Causality and the finite speed of wave propagation require that the 
amplitude must vanish before the direct arrivals of the transmitted component of the wave: 

K{x,t) = for t<x. (2.1) 

What we know about the wavefield is summarized in Table [T^. We now want to relate this 
equation to a focusing wave. A focusing wave is an incoming impulse wave incident on the 
medium from the left in our setup. Inside the medium it is followed by a 'tail' in such a way 
that at the focusing time t = the wave amplitude vanishes everywhere except at the focal 
point, where it gives rise to a 5-function peak. Such a cancellation of amplitudes rules out 
that refiected waves propagate to the left since otherwise they would keep propagating freely 
outside the medium even at focusing time. The wavefield of a focusing wave therefore has 
the form shown in Table [i|d. The unknown field in the medium is denoted by g'^{x,t). We 
next want to superimpose such waves so that they give rise to the wavefield of the refiection 
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problem outside the medium. This is done in Table [ij:;. For x ^ M [l^ and[T|:; are identical]^ 
Since the wave propagation solution is unique, the waves inside the medium must also be 



identical. Using eq. (2.1), we have: 



5{x 



t) + R{t + x)+ I g'{x, t') [6{t - t') + R{t + t')] dt' = for t < x. 

J — oo 



(2.2) 



which within the domain of validity is equivalent to, 

g%x, t) + R{t + x)+ [ g'{x, t')R{t + t')dt' = for t < x. 



(2.3) 



The reason why the upper bound of the integral is at a finite x deserves explanation: The 
bound has to be at least x so the waves are in agreement in x ^ M. It also has to be at 
most X to satisfy the constraint at t < x in x G M. The resulting equation (2.3) is the 



Gelfand-Levitan integral equation which allows solving for the focusing wavefield for a given 
reflection response R{t). 



2.2. Derivation in higher dimensions 

The derivation for higher dimensions proceeds along the same lines as in ID. The incident 
wave now enters through a surface S bounding the medium from one side as shown in Fig. [T] 
In principle the reflection problem can be generalized in different ways. The incoming wave 
can be an impulse wave incident at only one point on the surface S. This is a typical situation 
arising in quantum scattering and equations for a generalized Gelfand-Levitan equation have 
been derived for example in [12]. Particularly in seismic experiments, one deals with more 
than one point of incidence. Reflection responses are gathered all over a surface region S 
from shots all over that same region S. The reflection response -R(g, s, t) - the shot gather 
- therefore has two parameters: The position of the source s G S* and the receiver g G 5* 
in contrast to the one-dimensional case. This allows us to illuminate the medium from the 
entire surface S instead of from only one point. The incoming direct wave, which was a 
delta-function in ID, can therefore be generalized to an incoming wavefront passing through 
the entire surface S. In analogy to the ID wave which focuses at (x,t) = (0, 0) if it were to 
continue propagating freely, the Green's function G'^(x, —t) also focuses at the origin. The 
negative t is because it is the time reverse of the usual Green's function C^^Xjt), describing 
the propagation from a point source at the origin to x, which is the solution of, 

L-G'^(x,t) =5(x)(5(t), (2.4) 

where L is the operator of the wave equation. Note that we are not using its explicit form 
anywhere so results will be valid for all systems invariant under time reversal. When the wave 

^Strictly speaking, the product of delta distributions in the derivation is a dehcate issue, but in real- 
world situations the incoming wave cannot be a delta-distribution anyway. To derive this more rigorously, 
one would convolve the equations with an arbitrary source wavelet before superimposing. Then S{x — t') 
would be replaced with the convolved wavelet while 5{t — t') remains a delta function. Since the wavelet is 
arbitrary, it can afterwards be dropped again. 
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(a) Wavefield for the reflection problem: 

X ^ M 


X G M 


G"^(x, -t)+ / / i?(x, s, t - t')G'^{s, -t')dt'ds 

J SJ—oo 


K(x, t) with ir(x, t) = for t < t/ 


(b) Wavefleld of a focusing wave: 
X 


X e M 


G"'(x, -t) 


G^(x, -t) + G"(x, -t) 



(c) Wavefleld of the reflection problem in terms of focusing waves: 

X ^ M X G M 



G"'(s, -t') [(5(x - s)(5(f - t') + i?(x, s, < - t')] dt'ds 



S'/ — CXJ 



G"'(s, -t') [6{x - s)d{t - t') + i?(x, s, t - t')] dt'ds 



+ G''{s,-t')[S{x- s)S{t- t')+R{yi,s,t- t')]dt'ds 



G%x,-t)+ / / R{:>c,s,t-t') ■ G''{s,-t')dt'ds 

J S 'J —oo 



Table 2: In (c) we again superimpose focusing waves in such a way that the field is identical 
to (a) in the region x ^ M . 



is a time-reversed Green's function, the scattering field takes the form given in Table |2](a). 
The focusing wave is given in gb). And in|2{c) we again take a superposition of focusing 
waves so that they sum up to the wavefield of the scattering problem. In the ID case, the 
upper bound of the integral was x, assuming that the wave propagates with speed c = 1. 
We now have waves originating from all over S", and the integral bound is given by the first 
arrivals tj(x, s) from the source s to x. This is also a generalization from before, since we 
no longer have to assume constant speed of propagation. Apart from this difference, the 
argument proceeds as before and we compare [2|^a) and |2](c) in a; G M. For t < tj one 
obtains. 



G'^(x,-t)- 




R{:s^,s,t -t')G'^{s,-t')dt'ds+ / R{^,s,t - t')G'{s, -t)dt'ds = 0. 



SJ-oo 



(2.5) 



Since this equation is valid only for t < tf, the direct wave G"^(x, —t) does not contribute 
anything and has been dropped from the equation. This is in analogy to the dropped 6{t — x) 



in ID. Equation (2.5) is the 3D generalization of the Gelfand-Levitan equation. 



3. Proof of the iteration algorithm 
3.1. Integral equation from iteration 



We now derive the integral equation (2.5) from the iterative algorithm of Wapenaar et 
al. [7]. The algorithm relates downgoing waves p'^{g,t) and upgoing waves p~{g,t) through 
the following iteration, 

p+(g, -t) = G^g, t) - w{g, t)p-_,ig, t) , (3.1) 
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p~{g,t)= [ I R{g,s,t-t')p+{s,t')dt'ds 

J S J — oo 

The window function w(g,t) satisfies, 



w(g,t) 





for 


\t\ 




for 


\t\ 



(3.2) 



(3.3) 



At the point of convergence, we can insert eq. (3.1) into eq. (3.2) and obtain an integral 
equation for p~(g,t), 

I* poo r /»oo 

p~{g,t) =11 Rig,s,t-t')G'^{s,-t')dt'ds- / / R{g,s,t-t')w{s,-t')p'{s,-t')dt'ds 

JsJoo Js Joo , ^ 

foo p ptfls.r) 

R{g,s,t-t')G'^{s,-t')dt'ds- / / R{g,s,t-t')p-{s,-t')dt'ds 




S J ~oo 



S J ~tf {s,r) 



We can identify. 



G^(s,t) = -Po(s,t), 



(3.4) 
(3.5) 



and recover eq. (2.5) 



0=G'{g,t)+ / R{g,s,t-t')G'^{s,-t')dt'ds+ / R{g,s,t-t')G'{s,-t')dt'ds. (3.6) 



The integral bounds seemingly differ from eq. (2.5), but in the respective regions the contri- 



bution from the integral vanishes. The connection of the iteration algorithm with eq. (3.6) 
has already been proposed in [7] in eq. (8)-(9) and by [9j. 

3.2. Iteration from integral equation 



One way [15j to solve an equation of the form of eq. (2.5) is by introducing an auxiliary 
parameter A which we later set to unity: 



^^(g,^) = - f [ R{g,s,t-t')G\s,-t')dt'ds- X [ [ ^'"\{g,s,t-t')G'{s,-t')dt'ds 

JsJ Js J-tf {s,r) 



Writing G'^(g,t) as a power series, 

oo 

G\g,t) = J2m+iis,t), 

1=0 

we can solve for each order of A and obtain for A; G No, 

-.t^(s,r) 



(3.7) 




R{g,s,t - t')Gl{s, -t')dt'ds 



(3.8) 



S J ~tf (s,r) 
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where the starting value of the iteration is set to, 



G^(g,t) = G'^(g,t). (3.9) 

This is equivalent to the Wapenaar algorithm. The Banach fixed point theorem ensures 
that the integral equation has a unique solution and that the series converges for physically 
reasonable shot gathers. 

Before I conclude this section, let me briefiy remark on noise. The structure of the recursion 



eq. (3.8) as well as the iteration algorithm work is such that they provide an intuitive 
understanding of the noise to be expected in the data. Suppose that the refiection response 
-R(g, s, t) has been measured from an experiment where the incoming wave is a bandlimited 
source- wavelet s{t) - typically a Ricker function instead of an ideal 5-distribution. If we 
do not deconvolve -R(g, s), an error will be picked up at every iteration step. At every 
iteration, the new contribution would be convolved with the source wavelet s{t). The k-th 
contribution G^(g, t) is convolved /c-times with the source wavelet. In reality one should of 
course deconvolve R{g,s,t), but for a noisy signal deconvolution is never exact. A wave- 
front normally narrows down a bit after deconvolution but a sharp 6 peak will never be 
reached. Therefore qualitatively the argument remains true. In all cited papers on the 
subject, synthetic data has been used where deconvolution problems do not occur. Another 
way to think of noise is by realizing that the contribution at the k-th level comes roughly 
speaking comes from a wave which has propagated in the medium k times as long as the 
wave at iteration 1. 



4. A tour on scattering equations and their application 

A number of similar equations exist in the literature and have been applied to various 
scattering problems. It is worthwhile clearing up some confusion here by listing the main 
equations and their purpose. The one-dimensional Newton-Marchenko equation reads. 



g%e,x,t) 



/oo 
R{e, e', t' + e'x)g'{e', x, t')dt' (^4^) 

e-=-i,i e-=-i,i 



where g'^{e,x,t) = for t < ex. The parameter e = ±1 distinguishes between the direction 
of wave propagation. Rose initially used this equation in his work and showed that it solves 
for the field g^{e,x,t) that focuses to 6{x) at t = Its 3D generalization is given in 
slightly different form in [11] and can be written, 

g'{ei,tf,t)= [ R{ei,es,t + tf)des+ [ [ R{e,,e' ,t + t')g%-e' ,tf,t')de'dt' . (4.2) 

where g^{ei,tf,t) is the scattered field for an impulse injected in direction ej. In eq. (4.1) 
the refiection response is recorded on both sides of the medium, e' = 1 and e' = — 1. In 
other words, the equation depends on the refiected as well as the transmitted wave. In 
eq. (4.2) the summing over both sides has been replaced by an integral over a closed surface 
S"^ surrounding the entire medium. The scattered wavefield is recorded all over the closed 
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Figure 1: A schematic experimental setup with sources and receivers on the surface S. 
With the proposed method (L) one-sided surface measurements are sufficient. Seismic 
interferometry (m.) needs a source inside the medium for each subsurface point of interest. 
The Newton-Marchenko equation requires receivers surrounding the entire medium. 



surface. These equations are known from scattering theory, where they have been apphed 
to the plasma equation. In momentum space it reads, 



{A + k'^ -V{x))u{k,x) =0. 



(4.3) 



Note also that in the scattering problem an impulse signal (or in momentum space a plane 
wave) is scattered as it passes through a medium and the scattering response is recorded. This 
wave is incident from one specific direction only, here e^. From the theoretical point of view 
it is of course attractive that a simple plane wave is all that is needed to solve the scattering 
problem, but in the presence of noise it is preferable to have an incident wavefield coming from 
an entire area covered with sources. Furthermore, for seismic experiments it is impractical to 
bury sources or receivers in the ground. We want to deduce the properties of the subsurface 
from measurements on the surface only. In other words, we are looking for an equation 
using one-sided illumination only. Rose has solved that problem in ID using focusing waves. 
In ID the incident wave would be a delta function impulse followed by a tail such that 
all scattered waves are cancelled at the focusing time. It was shown that an equation 
called Marchenko equation or Gelfand-Levitan equation solves the problem 03] • While 
for example in [10] the authors refer to Newton and Marchenko by dubbing their imaging 
method Newton-Marchenko-Rosen imaging, their single-sided illumination algorithm in ID 
actually corresponds to the Gel'fand-Levitan equation. It has been rederived in this paper 



and is given in eq. (2.3). The Newton-Marchenko and the Gel'fand-Levitan equations might 



look similar but they are rather different. The former needs illumination from all sides, 
the time integral is unbounded and the solution is the entire wave-field. The latter needs 
illumination only from one side, the integration ends at the time of the first arrivals and 
its solution is only the scattered field. The equations appear in the literature in somewhat 
varying forms and also for different applications. Characteristic of Marchenko-type equations 
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Common name 



Eqn Dim recording surface 



Newton-Marchenko 
3D Newton-Marchenko 
Gelfand-Levitan or Marchenko 
N/A 



(4.1) 



(4.2) 



(2.3) 



(2.5) 



ID both sides 

3D closed surface 

ID single-sided 

3D single-sided 



are integral limits of only ±oo or in the time integration whereas Gelfand-Levitan equations 
have one integration bound at a different value. A " Generalized Gelfand-Levitan" equation 
also exists in the literature [TT]. Like all Gelfand-Levitan type equations it has a finite 
bound of integration over the time variable, but it also requires an integration over a closed 
surface around the medium, so it is not a single-sided scattering equation. Its integration 
kernel is also different from the ones in the table. To avoid causing confusion, I therefore 
abstained from calling the equation in Table |4] generalized Gelfand-Levitan equation or 3D 
Gelfand-Levitan equation. 



5. Summary and conclusion 

I have derived an integral equation equivalent to the Wapenaar iteration algorithm from 
simple physical reasoning only. The equivalence with the iteration algorithm was proven. 
The formula can be applied for different purposes. One can find a Green's function for a 
'virtual source' under the subsurface, which is nothing but the time-reversal of the wavefield 
of the focusing wave ioi t > tj [H|7] . The advantage is that it has been recovered without 
the need to place a source or receiver inside the medium. From the fundamental identity 
(dubbed by Newton as "miracle equation") [TT] the potential V^(x) of the plasma wave 



equation (4.3) can be recovered. And by separating the up-going from the down-going wave, 
the Green's function can be used as the input of an imaging condition and provide an image 
with all internal multiples. Examples with synthetic data are shown in [TU]. Using the 
integral equation, instead of the Wapenaar algorithm there are presumably more effective 
ways to compute images, perhaps using a layer-peeling algorithm. A variety of applications 
are thinkable since knowing the Green's function for every point inside a medium amounts 
to knowing everything about the scattering problem. 
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